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Bilayer graphene (BLG) with a tunable bandgap appears interesting as an alternative to graphene 
for practical applications, thus its transport properties are being actively pursued. Using density 
functional theory and perturbation analysis, we investigated, under an external electric held, the 
electronic properties of BLGs in various stackings relevant to recently observed complex structures. 

We established the hrst phase diagram summarizing the stacking-dependent gap openings of BLGs 
for a given held. We further identihed high-density midgap states, localized on grain boundaries, 
even under a strong held, which can considerably reduce overall transport gap. 

PACS numbers: 61.48.Gh, 73.22.Pr, 73.21.Ac 


The discovery of graphene has opened new avenues for 
studying the role of dimensionality on the fundamental 
properties of materials [T]. Although graphene shows 
excellent electrical properties [2], the zero bandgap of 
graphene limits its practical application as an electronic 
device. On the other hand, gap opening is possible in 
BLG, thus making it a very promising material that over¬ 
comes graphene’s key limitation while retaining many of 
its interesting properties. For example, massive Dirac 
fermions in BLG exhibit a bandgap tunable by apply¬ 
ing a transverse electric held (E-held) [3]; this has been 
demonstrated by optical [4] and electrical transport mea¬ 
surements using dual-gated devices 0IH1- However, these 
measurements leave a couple of unsolved problems: 1) 
the origin of unexpectedly small transport gaps that are 
two orders of magnitude smaller than optical gaps [5] 
and 2) the origin of anomalous low-temperature (< 2 
K) transport behaviors dominated by hopping between 
localized midgap states, presumably induced from disor¬ 
ders or defects urn- 

Recent experiments have revealed complex conhgura- 
tions in BLG, including various stacking domains induced 
by rotational faults and soliton formation jSUTO]. While 
AB stacking is energetically most favorable, the non-AB- 
stacking region can be stabilized by a minute twist m 
and the stacking boundary [8] . The local stacking conhg- 
uration is strongly coupled to its electronic structure and 
its response to an external E-held. Therefore, it is criti¬ 
cally important, fundamentally and practically, to under¬ 
stand the observed complex stackings and their impact 
on the overall electronic properties. 

In this letter, using the framework of an effective 
Hamiltonian based on density functional theory (DET) 
and perturbation theory, we analyze gap-opening prop¬ 
erties of BLGs near the high-symmetry stackings (AA, 
AA', and AB), under an applied E-held. We establish a 
phase diagram for the stacking-dependent gap openings, 
and further identify grain boundaries containing non-AB 
stackings as a source for high-density midgap states even 


under a strong E-held. Our hndings offer insight to un¬ 
derstanding the intrinsic transport properties of BLGs. 

Our DET calculations adopt the Perdew-Burke- 
Ernzerhof version of exchange-correlation functional m 
and the projector augmented wave method m for ionic 
potentials as implemented in the Vienna Ab Initio Sim¬ 
ulation Package M- We obtain interlayer distances be¬ 
tween 3.25 A (AB) and 3.45 A (AA) with van der Waals 
correction m-, interlayer distance of all the conhgura- 
tions is hxed at 3.35 A (unless specihed) with practically 
no changes in their band structures. To ensure an ac¬ 
curate bandgap, the 2D DET band structure near the 
K point is interpolated |16] using maximally localized 
Wannier functions m- Effective Hamiltonians are con¬ 
structed with the obtained hopping parameters truncated 
to the hrst nearest interlayer hoppings (see details in Sup¬ 
plemental Materials). 

One of the intriguing properties of BLG is that a 
change in weak interlayer interaction (which is an order of 
magnitude smaller than the intralayer coupling strength) 
accompanied by a modihcation in stacking conhguration 
can signihcantly alter the electronic structure around the 
Eermi level. Eigure[^a) illustrates schematic band struc¬ 
tures of AA, AA', and AB stacking, where we dehne sys¬ 
tems with equivalent two sublattices, such as AA and 
AA', as sublattice-symmetric systems; otherwise, they 
are sublattice-asymmetric stackings (such as AB). 

Eigure [^b) shows the atomistic modeling of an ex¬ 
perimentally observed domain boundary [9]. Here, for 
visual clarity, we made the phase boundary region (the 
region of continuous structural transition between two 
AB-stacking regions) much smaller than the experimen¬ 
tal ones. Eigure [^c) plots stacking-dependent poten¬ 
tial energy with optimized interlayer distances in the 2D 
translation vector space, where AB stacking is used as 
a reference point. The arrow denotes the displacement 
vector between the left and right domains in Eig. Bb)- 
Local stacking configurations of the transition region are 
distributed on this arrow. To remove this soliton-like 
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FIG. 1. (Color online), (a) Schematic band structures of AA, 
AA', and AB. The solid lines are reflection planes, where the 
translation vectors [Vx^Vy) describe the relative displacement 
between the two layers in the xy plane, (b) Modeling of the 
two AB-stacking boundary, (c) Stacking-dependent potential 
energy of BLG per unit cell, where the origin corresponds to 
AA stacking. A lattice Wigner-Seitz cell is highlighted by 
the solid white line, and the arrow denotes the displacement 
vector between the two AB-stacking domains shown in (b). 
(d) Minimum energy path between the two AB stackings. 


boundary, one needs to displace the one on the left or 
right domain by a displacement vector. The minimum 
energy path between the two AB stackings lies along 
the edge of the hexagon with an energy barrier of 5.3 
meV/cell (see Fig. [^d)). Though this energy barrier 
seems quite small, the stacking domain should move as 
a whole so that the high energy barrier, which is pro¬ 
portional to the area of the domain (> 10^ unit cell), 
should be overcome. This explains the observed stability 
of non-AB stacking regions. 

The gap-opening mechanism of BLGs can be highly 
stacking dependent. Thus, we first examine the individ¬ 
ual band structures near the high-symmetry stackings. 
We then discuss the gap-opening properties across com¬ 
plex domain boundaries. Their effective Hamiltonian in 
crystal momentum {k) space can be described by a 4 x 
4 matrix with the basis Aup^ Bup, Adn^ and where 
A and B represent sublattice indices and the subscripts 
up and dn denote the upper and lower graphene layer, 
respectively. Here, the 2x2 block-diagonal components 
correspond to the individual graphene layers while all 
others describe interlayer coupling. We will now focus 
only on the effective Hamiltonian near the K points; 
band structures around K' can be obtained by applying 
time-reversal symmetry to those of K. 

First we consider the configurations of A A stack¬ 
ing. Neglecting small Bloch phase variations under the 
Fourier transformations of interlayer coupling, we find 


that the Hamiltonian of AA stacking around K becomes 


H^{k) + 


0 0 7aa 0 

000 7aa 

lAA 000 

0 7aa 0 0 


with Ho{k) defined as 


Hoik) = 


0 hupk-^ 0 

hvpk- 0 0 

0 0 0 

0 0 hvpk- 


0 

0 

hupk-^ ’ 

0 


( 1 ) 


( 2 ) 


where the Fermi velocity multiplied by the reduced 
Planck constant becomes hvp = ^ ^ 5.4 eV-A and 
k± = ky Al ikx. JAA {= —0.34 eV) is obtained by the 
Fourier transformation of the interlayer hopping between 
Aup (Bup) and Adn (Bdn), Iaa [H] • The hopping 
parameters between Aup {B^p) and B^n (Adn) become 
zero because the Bloch phases of three interlayer near¬ 
est neighbors cancel each other at the K point; that is, 
lAB = 0. 

By changing our basis to the bonding and antibonding 
state of each sublattice, the decoupling of two Dirac cones 
becomes more transparent as follows: 


Hoik) A 


lAA 0 0 0 

0 7aa 0 0 

0 0 -7aa 0 ’ 

0 0 0 -7aa_ 


( 3 ) 


which is a block-diagonal Hamiltonian describing two 
Dirac cones with energy shift ±7 aa as shown in the 
schematic band diagram of Fig. [^a). 

In the AA'-stacking configuration, one can also ex¬ 
plicitly illustrate the decoupling of Dirac cones by 
changing the basis to interlayer bonding and anti¬ 
bonding of phase-shifted sublattices [^AupiBup) A 

exp i^^)AdniBdn)]- The Hamiltonian of A A' stack¬ 
ing then becomes 


Hoik) A 


lAA —lAB 0 0 

-"fAB lAA 0 0 

0 0 -^AA lAB 

0 0 7ab -7aa 


( 4 ) 


where ^aa — —0.11 eV and ^ab — —0.22 eV, corre¬ 
sponding to two Dirac cones separated by 0.22 eV in 
energy with an additional 0.08 splitting in /c-space. 
Wavefunctions of the decoupled Dirac cones of both AA 
and AA' stackings have interlayer antibonding and bond¬ 
ing characteristics, depicted respectively in red (shaded) 
and blue (hatched) in Fig. [^a). 

The Hamiltonian of AB staking can be written as 


Hoik) A 


0 00 ^ab 

0 0 0 0 
0 0 0 0 ’ 
^ab 00 0 


( 5 ) 
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FIG. 2. (Color online) Projected band structures around K point. Each configuration is represented by the translation vector 
in the irreducible zone of the lattice Wigner-Seitz cell (triangle), where the lower vertex dehnes the origin. and k± are 
dehned for each conhguration. Energy (/c-space) ranges from —0.5 eV (—0.2 A“^) to 0.5 eV (0.2 A“^) relative to the Eermi 
level (K). Band structures near the high-symmetry stackings are projected onto the /cy-energy and /c^-energy planes without 
and with an E-held. The inset in the third column of AA stacking highlights a small bandgap (~ 10 meV). 


that is, there are doubly degenerate states at the Fermi 
level composed of one sublattice per layer, with no direct 
coupling between them. In this case two Dirac points are 
merged at the K point and split into bonding, antibond¬ 
ing, and nonbonding types (see Fig. [^a)). 

Next, we trace how a small translation or an external 
E-field can change the band structures, especially near 
the Fermi level. Our approach is to treat small trans¬ 
lations and E-fields as perturbations to the individual 
high-symmetry stackings. 

The A A panel of Fig. summarizes projected (onto 
k\\ and k±) band structures for around-AA-stacked 
graphenes without and with an E-field. Since interlayer 
hopping parameters {jaa and ^bb) are the same, one 
cannot generate an onsite energy difference in the 2x2 
diagonal block simply by atomic translation, which ex¬ 
cludes the direct coupling between two crossing bands 
(z.e., no gap opening). At the Fermi level, the hole band 
of one Dirac cone is degenerated along with the electron 
band of the other Dirac cone. The Fermi surface of the 
AA stacking is the intersection of two vertically (energet¬ 
ically) shifted cones, a circle. A small translation results 
in a slight /c-shift and a coupling of two Dirac cones; a k- 
shift changes the circular intersection into a tilted ellipse 
(which is still a circle when projected on /c-space), while a 
coupling introduces energy splitting at the intersection. 
In general, the energy splitting depends on the angu¬ 
lar position of the intersection and becomes zero at two 
points. These form two crossing points near the Fermi 
level as shown in the second row. An applied E-field in¬ 
troduces an additional energy splitting that also depends 


on the angular position and becomes zero at two points. 
When an E-field is combined with a sublattice-symmetric 
translation, their zero splitting points coincide and the 
system remains metallic. In contrast, with the sublattice- 
asymmetric translation, each zero coupling points are at 
a different position and the crossing points disappear as 
shown in the inset in the fourth row. Especially, when 
sublattice-asymmetric translation is applied toward AB 
stacking in the presence of a reflection and time-reversal 
symmetry, the minimum bandgaps occur along k\\ and 
are located exactly at the same energy. This means that 
the critical field for opening a gap is infinitesimally small. 
The perturbational results on the size of the bandgap are 
summarized in Table 1. As an example, the fourth row 
in the A A block of Fig. shows a small bandgap of ^10 
meV (see the inset) for an asymmetric translation of 0.3 
A and an E-field of 0.5 eV/A. 

Changes in band structures for around-AA' stacking 
are well pronounced in the k\\ =0 plane [the blues lines 
in the second and fourth rows of the A A' block in Fig. 
1^. Of the four bands in that plane, only different Dirac 
cones can be coupled by a translation. In contrast, un¬ 
der sublattice-symmetric translation, only parallel-band 
pairs of each Dirac cone are coupled, resulting in a bal¬ 
anced repulsion between them. On the other hand, un¬ 
der sublattice-asymmetric translation only non-parallel- 
band pairs of each cone are coupled, which induces an 
unfavorable crossing. Also in this slice, E-field only cou¬ 
ples parallel-band pairs for sublattice-symmetric transla¬ 
tion. Though the crossing point in a Dirac cone does not 
open, each Dirac cone’s crossing band now has a small 
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TABLE L Analytic expressions of (pseudo-) gaps when a small translation x from reference stacking configurations is combined 
with an interlayer potential difference f/, where ADE(k) denotes energy (crystal momentum) separation of two Dirac points. 


reference 

AA 

AA' 

AB 

translation 

direction 

toward AB 

toward AB 

any direction 

(pseudo-) gap 

. ADk{x) jj 

2ADEix)^ 

A7(x) 

ADe{x)'^ 

^AB jt 


parameters 

hvp = 5.4 eV-A 

ADk{x) = 0.03a: 
ADe{x) = 0.68 eV 

Aa{x) = Re ^exp(—^z) 

0.2x < A 7 (x) < 0.3x (eV/A) 
ADe(x) = 0.22 eV 

^AB — 0.30 eV 

U {E = 0.5 V/A) 

0.15 eV 

0.52 eV 

0.55 eV 


component of the opposite Dirac cone. If we apply an E- 
field with sublattice-symmetric stacking, crossing bands 
still remain crossed because in a Dirac cone, one crossing 
band does not have a component parallel to the other 
crossing band. But if an E-field is applied to sublattice- 
asymmetric stacking, each crossing band now has a small 
component parallel to the other crossing bands, which 
opens a small bandgap. In spite of Dirac cones open¬ 
ing, the energy level of each Dirac point is different [the 
fourth row in the AA' block of Fig. [^, thus a relatively 
strong E-field is required to change this pseudogap into 
a true gap. 

Finally, we move on to the properties of around-AB 
stacking. As the stacking deviates from exact AB, the 
doubly degenerate states of AB at the K point imme¬ 
diately split into two crossing points. This feature is 
demonstrated in the second and third columns of the 
AB block in Fig. From a symmetry viewpoint, the 
threefold rotational symmetry of monolayer graphene is 
recovered in A A- and AB-stacked BLG. Combined with 
translational symmetry, this imposes a threefold symme¬ 
try around the K point. Because two separated crossing 
points are not compatible with the symmetry, wavefunc- 
tion symmetries change during the merging of two cross¬ 
ing points. It was reported that this merging process is 
very sensitive to miniscule translation (^0. 01 A) and the 
band topology near the Fermi level changes m- Around 
AB stacking, an E-field opens a bandgap. Especially, 
from the eigenvalues of the Hamiltonian, the bandgap is 
— 7 ab where ^ab = = 0.30 eV. All the 

perturbational results for the bandgap are summarized 
in Table 1. 

Figure |^a) presents the stacking-dependent bandgap 
under a perpendicular E-field of 0.5 V/A. A sizable 
bandgap opens only around the AB stacking while the 
rest still remains metallic. As E-field goes to zero, the 
metal-semiconductor phase boundary approaches the line 
connecting AA and AB stacking, and the entire region 
becomes metallic as shown in Fig. [^b). Though no 
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FIG. 3. (Color online) (a) A stacking-dependent bandgap 
under a perpendicular E-field of 0.5 V/A. A lattice Wigner- 
Seitz cell is shown by broken lines and an irreducible zone 
by solid lines, (b) Metal-semiconductor phase boundaries for 
different electric field strengths are shown for the irreducible 
zone, (c) Local stacking configurations of simulated structure 
are represented by colors in the triangle at left, (d) Local 
densities of states of the spotted region in (c) are plotted 
from the left with (dotted line) and without (solid line) an 
E-field. 


bandgap opens by a pure translation, a minute bandgap 
(<7 meV) was reported [20] for a specific rotation angle 
without any E-field. To investigate the effect of the non- 
AB stacking region on the transport property, we con¬ 
structed an atomic model of a stacking domain bound¬ 
ary with a transition length of 50 A (Fig. § (c)). Tight- 
binding parameters are assigned to each atom according 
to its local stacking configuration [21]. When 0.5 eV of 
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onsite energy difference between two layers is applied, 
which corresponds to 0.5 V/ A of E-field, a bandgap opens 
at the AB stacking region while there remains finite den¬ 
sity of states at the non-AB stacking region (Fig. |^d)). 
This indicates that a high density of midgap states is 
localized along stacking boundaries even under a strong 
E-field. Because the apparent transport gap is actually 
estimated from the activation energy of the carrier, a 
conduction through these midgap states can explain the 
small transport gap and the low-temperature hopping 
transport in dual-gated devices. 

In summary, we theoretically investigated stacking- 
dependent gap-opening properties of symmetry-broken 
bilayer graphenes, and established a bandgap phase di¬ 
agram. Our findings may prove to be instrumental in 
developing graphene-based electronic devices. 


ACKNOWLEDGMENTS 

This research was supported by the Center for 
Nanophase Materials Sciences, Oak Ridge National Lab¬ 
oratory by the Scientific User Facilities Division, Of¬ 
fice of Basic Energy Sciences, U.S. Department of En¬ 
ergy. J. R., S. H., and G. K. were supported by 

(1) the NanoMaterial Technology Development Pro¬ 
gram (2012M3A7B4049888) through the NRF, funded 
by the Ministry of Science, ICT and Future Plan¬ 
ning; (2) the Priority Research Center Program (2010- 
0020207); and (3) the Basic Science Research Program 
(2013R1A2009131) through NRF, funded by the Min¬ 
istry of Education in Korea. 


Corresponding Author: gunnkim@seiong.ac.kr 
^ Corresponding Author: niyoon@ornl.gov; Notice: This 
manuscript has been authored by UT-Battelle, LLC, un¬ 
der Contract No. DE-AC05000R22725 with the U.S. De¬ 
partment of Energy. The United States Government re¬ 
tains and the publisher, by accepting the article for publi¬ 
cation, acknowledges that the United States Government 
retains a non-exclusive, paid-up, irrevocable, world-wide 
license to publish or reproduce the published form of this 
manuscript, or allow others to do so, for the United States 
Government purposes. 

[1] K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. 
V. Khotkevich, S. V. Morozov, and A. K. Geim,, Proc. 
Nat. Acad. Sci. 102 , 10451 (2005). 

[2] K. I. Bolotin, K. J. Sikes, Z. Jiang, M. Klima, G. Fuden- 
berg, J. Hone, P. Kim, and H. L. Stormer, Solid State 
Commun. 146, 351 (2008). 

[3] H. K. Min, B. Sahu, S. K. Banerjee, and A. H. MacDon¬ 
ald, Phys. Rev. B 75, 155115 (2007); E. McCann, Phys. 
Rev. B 74, 161403(R) (2006). 

[4] T. Ohta, A. Bostwick, T. Seyller, K. Horn, and E. Roten- 
berg. Science 313, 951 (2006); K. F. Mak, C. H. Lui, J. 
Shan, and T. F. Heinz, Phys. Rev. Lett. 102, 256405 


(2009); Y. B. Zhang, T. T. Tang, C. Girit, Z. Hao, M. 

C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen, and F. 
Wang, Nature (London) 459, 820 (2009). 

[5] T. Taychatanapat and P. Jarillo-Herrero, Phys. Rev. 
Lett. 105, 166601 (2010). 

[ 6 ] J. B. Oostinga, H. B. Heersche, X. Liu, A. F. Morpurgo, 
and L. M. Vandersypen, Nature Mater. 7, 151 (2008); T. 
Taychatanapat and P. Jarillo-Herrero, Phys. Rev. Lett. 
105, 166601 (2010); J. Velasco Jr, Y. Lee, Z. Zhao, L. 
Jing, P. Kratz, M. Bockrath, and C. N. Lau, Nature Nan- 
otechnol. 7, 156 (2012); A. Varlet, M. H. Liu, V. Krueckl, 

D. Bischoff, P. Simonet, K. Watanabe, T. Taniguchi, K. 
Richter, K. Ensslin, and T. Ihn, Phys. Rev. Lett. 113, 
116601 (2014). 

[7] S. Tanabe, Y. Sekine, H. Kageshima, M. Nagase, and H. 
Hibino, Jpn. J. Appl. Phys. 50, 04DN04 (2011). 

[ 8 ] J. Lin, W. Fang, W. Zhou, A. R. Lupini, J. C. Idrobo, J. 
Kong, S. J. Pennycook, and S. T. Pantelides, Nano Lett. 
13, 3262 (2013). 

[9] J. S. Alden, A. W. Tsen, P. Y. Huang, R. Hovden, L. 
Brown, J. Park, D. A. Muller, and P. L. McEuen, Proc. 
Nat. Acad. Sci. USA 110 , 11256 (2013). 

[10] H. Hibino, S. Mizuno, H. Kageshima, M. Nagase, and H. 
Yamaguchi, Phys. Rev. B 80, 085406 (2009). 

[11] K. S. Kim, A. L. Walter, L. Moreschini, T. Seyller, K. 
Horn, E. Rotenberg, and A. Bostwick, Nature Mater. 12 , 
887 (2013). 

[12] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. 
Lett. 77, 3865 (1996). 

[13] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 

[14] G. Kresse and J. Furthmuller, Phys. Rev. B 54, 11169 
(1996). 

[15] A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102 , 
073005 (2009). 

[16] We use extremely dense grid points (0.001 grid spac¬ 
ing) to ensure an energy resolution up to 5 meV. Max¬ 
imally localized Wannier functions are constructed from 
the Bloch functions of uniformly-sampled 36 x 36 x 1 
k points in the Brillouin zone. In the energy range of 
±2 eV from the Fermi level, first-principles band struc¬ 
tures at the sampled /c-points are exactly reproduced in 
this interpolation scheme. 

[17] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 
(1997); 1. Souza, N. Marzari, and D. Vanderbilt, Phys. 
Rev. B 65, 035109 (2002). 

[18] A real space interlayer hopping parameter 7 is related to 
/c-space one 7 by 

~ V ^ i-fC-(R- ttt,—R-7 ) 

= 2^ e ' ™ ‘'7(!,R,)(m,R™)> 

where /, m denotes sublattice indices and Rz(m) de¬ 
notes position vector of l{m) sublattice atom and K = 
(|^, ^^) and a is G-G bond length. 

[19] Y.-W. Son, S.-M. Ghoi, Y. P. Hong, S. Woo, and S.-H. 
Jhi, Phys. Rev. B 84, 155410 (2011). 

[20] S. Shallcross, S. Sharma, and O. A. Pankratov, Phys. 
Rev. Lett. 101, 056803 (2008). 

[21] Based on the uniformly sampled 487 stacking configura¬ 
tions, we extract the exact hopping parameters and trun¬ 
cated them up to nearest-interlayer pairs. The grid data 
of hopping parameters is used in the linear interpolative 
mapping to each atom with the same local stacking con- 



6 


figuration. Two hundred k-points along the boundary are 


sampled for the density of states calculations. 



